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Abstract 

The onset of convection in a horizontal layer of fluid heated from below in the 
presence of a gravity field varying across the layer is investigated. The eigenvalue 
problem governing the linear stability of the mechanical equilibria of the fluid layer 
in the case of free boundaries is solved using a Galerkin method based on shifted 
polynomials (Legendre and Chebyshev polynomials). 
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1 Statement of the problem 

Physical problems concerning the motion of fluids in the presence of a variable grav- 
ity field can be encountered in many practical applications, e.g. convection problems 
in porous media, crystal growth domain or problems concerning the mass transport in 
the Earth's system. That is way, convection problems with variable gravity fields have 
been intensively studied. In this paper we analyze the influence of such gravity fields 
varying across the layer on the stability bounds in a convection problem that arises in a 
horizontal layer of fluid heated from below. The gravity field acting in the ^-direction is 
orthogonal to the fluid layer and is assumed to depend on the vertical coordinate z only 
[12j . For such a variable gravity field different points of the fluid experience different 
buoyancy forces. As a consequence, part of a fluid layer tends to become unstable while 
the other tends to remain stable, the mechanical equilibrium turning into a convective 
motion. 

Experimental measurements of the Earth's upper atmosphere show that the atmo- 
spheric density decreases as the altitude increases as an approximately exponentially 
function of the vertical height. 

Taking these into account, consider a layer of heat-conducting viscous fluid con- 
tained between the planes z = and z = h[l2\. The equations governing the convective 
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motion and the conducting state are [8] 



<9v i 

— + (v • grad)v = — -gradp + i/Av + g(z)aT, 

divv = 0, t > (1) 

dT 

— + (v • grad)T = kAT, 

where v is the coefficient of kinematic viscosity, p the density, a the thermal expansion 
coefficient, k the thermal diffusitivity, p is the pressure, T the temperature, v the 
velocity and the gravity g{z) is defined by g(z) = gH(z)\t, where g is a constant. 

The linear stability of the conduction stationary solution characterized by v = of 
equations (1), written in the nondimensional form, against normal mode perturbations, 
is governed by a two-point problem for the ordinary differential equations 

(D 2 - a 2 ) 2 W = RH(z)a 2 e, 

(D 2 - a 2 )@ = -RN{z)W, { ' 

where D = — , R 2 is the Rayleigh number, a is the wavenumber and W and G are the 

dz 

factors in v and 6 respectively, depending on z. 

Consider N(z) = 1 and H{z) = 1 + eh(z), z £ (0, 1). The parameter e represents 
a scale for h{z). In this case, the two-point problem for (2) consists of the ordinary 
differential equations [12] 



(D 2 - a 2 ) 2 W = R{1 + eh(z)]a 2 e, 
(D 2 - a 2 )G = -RW 



(3) 



and the usual boundary conditions for free boundaries read 

W = D 2 W = B = at * = 0,1. (4) 

We look for the smallest eigenvalue R (the Rayleigh number) in ©-(HI) defining the 
neutral manifold. 

In the case of rigid boundaries, i.e. the boundary conditions are 

W = DW = G = at z = 0, 1 (5) 

the eigenvalue problem has been studied by us in some previous papers, e.g. [3], [3]. 
The problem (E]) - <[5j) was also investigated in |12| . Straughan performed numerical 
evaluations of the Rayleigh number by using the energy method for some various gravity 
fields and our numerical evaluations obtained in [3], [JJ were similar to those obtained 
in [12] . Herein, for the same varying gravity fields we will obtain numerical results in 
the case of free boundaries. 

In [6], [7] it is proved that when H(z)N{z) > across the layer the principle 
of exchange of stabilities holds. In our case N{z) = 1, e is a small parameter, for 
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the chosen functions h{z) the inequality H(z) > is valid, therefore the principle of 
exchange of stabilities holds no matter what the boundary conditions. We can mention 
that there are cases in which this condition is not satisfied and the principle of exchange 
of stabilities holds. 

Although their study started in 1961, the eigenvalue problems governed by systems 
of ordinary differential equations with variable coefficients and depending on many 
physical parameters are difficult to solve. In most cases, an approximative solution 
for this type of problems can be obtained using spectral methods. Since it has the 
advantage of optimal analysis, i.e. the analytical study is not a tedious one and the 
obtained numerical results are very good compared to other methods, the Galerkin 
method was chosen for the study of the eigenvalue problem ©-Q. 

The expansion sets of functions used for the various fields encountered in the con- 
vection problems from hydrodynamic stability theory (e.g. the velocity field, the tem- 
perature field, the concentration field) must have a basic property: they must be easy 
to evaluate. That is why, most of the times, the trigonometric and the polynomial 
functions, easy to evaluate, are used. A second property is the completeness of the 
expansion sets of functions. This assures that each function of the given space can be 
written as a linear combination of functions from the considered set (or, more likely, 
as a limit of such a linear combination). The Chebyshev polynomials, the Legendre 
polynomials, the Hermite functions, the sine and cosine functions, satisfy this condition. 

In [3] the two-point problem ©-([I]) was investigated using methods based on Fourier 
series of trigonometric functions. Herein, the trial and the test sets of functions will 
consists in Chebyshev and Legendre polynomials. However, when the boundary condi- 
tions are very complicated the Galerkin approach is not easy to apply. That is way, in 
order to write the boundary conditions in a simpler form, we introduced the function 
^ = {D 2 - a 2 )W and taking into account that W = D 2 W = at z = 0, 1 we get ^ = 
at z = 0, 1. 

By denoting U = (W,^,Q) the eigenvector in ©-([I]), the two-point problem can 
be rewritten 

L{U = (D 2 - a 2 )W - * = 0, 

L 2 U = {D 2 - a 2 )^ - R(l + eh{z))a 2 e = 0, (6) 
L 3 U = 44> {D 2 — a 2 )Q + RW = 0, 

with the boundary conditions 

W = ¥ = e = 0at* = 0,l. (7) 

2 Methods based on shifted polynomials 

In order to perform not only an analytical study, but also to obtain numerical evalu- 
ations for the Rayleigh number R 2 , each of the functions from the unknown eigenvector 
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U, is approximated by a truncated series of orthogonal polynomials, in our case Cheby- 
shev and Legendre polynomials. These polynomials are orthogonal on [—1,1]. Since 
in this problem the range is [0, 1], we will use shifted polynomials, orthogonal on [0, 1], 
obtained from the original polynomials by a variable transformation. 

The Chebyshev polynomials were widely used in spectral methods for ordinary 
differential equations, e.g. pQ, [5], [8]. Here we present only some basic properties of 
these polynomials necessary for our study. 

The Chebyshev polynomials (of the first kind) of degree n, T n (z), are orthogonal on 

1 f 1 

[— 1, 1] with respect to the weight function w(z) = , i.e. / T n (z)T m (z)w(z) = 

VI — z 2 J-i 

7r r 2 if n = 

7rC n (5 mn , c n = < ' ' . The shifted Chebyshev polynomials (of the first kind) (SCP) 

of degree n on (0, 1), T*(z), are defined by the relation T*(z) = T n (2z — 1). The fol- 
lowing orthogonality relation holds 

1 T:(z)T^(z)w*(z)dz = ( 2 C "/™7 if 1 = j > (8) 
o [ 0, if % + j, 

with respect to the weight function w*(z) = — -. The recurrence relation between 

z{l-z) 

T* has the form 

T* n {z) = 2(2z - l)T* n ^(z) - T:_ 2 (z). 

Similarly with Shen[llj, let us introduce M\ = {§* k {z)}k&, a the complete set of 
orthogonal functions in L 2 (0, 1), $%(z) defined by 

<J?£(z)=T fc *(z)-T fc * +2 (z) (9) 

and satisfying boundary conditions of the type <&£(0) = 3>£(1) = 0. Then the un- 
known functions W, G can be expanded upon the complete set M\ and they satisfy 
automatically all the boundary conditions. We have 

n n n 

w=Y,w k n(z), e = ^e fc ^,(z), *=Y,*kn(z). m 

k=0 k=0 k=0 

The system ([3]) can then be written in terms of the expansion functions only. Im- 
posing the condition that left-hand side equations of the system to be orthogonal on 
£ = 0,1, ...,n, we get the algebraic system 



' t {{{D 2 - a 2 )^{z)^(z))w k - {^t(z)^(z))^ k } =^ 

± { ({D 2 - a 2 )$* (z), $*(z)) * fe - i?a 2 ((l + eh(z))^ k (z)^*(z))e k \ = 0, (11) 
E { ((D 2 - a 2 )$l(z), Q k + R(& h (z), **(z))W k } = 



k=0 
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in the unknown coefficients Wk, VPfc, Ofc- Since not all these coefficients are null, the 
condition that the determinant of the system vanish is imposed leading to the secular 
equation. 

Following [8] it is easy to deduce the following derivation formulae 

fc-i fe+i 
(*;(*))' = 2{2fc T r *(z)-2(k + 2) T r( z )} ( 12 ) 

r = r = 

k — r odd k + 2 — r odd 

for the first derivative of the function and 

fc-2 fc 
(*£(*))" = 4{ ^ (fc-r)fc(*;+r)^(2!)- (A;+2-r)(A:+2)(A:+2+r)T;(z)} 

r = r = 

k — r even fc + 2 — r even 

(13) 

for the second one. In the numerical evaluations we will take into account that the first 
term in each of the involved sums is halved. 

The presence of the varying gravity field led to nonconstant coefficients, such that 
the analytical expression of the scalar product (h{z)<&* k {z) , <&*(z)) is based on the rela- 
tion [5] 

z r T s {z) = ±J2ciT s - r+2i (z). (14) 

i=0 

The analytical expressions of all the other scalar products from (jlip are taking with 
respect to the weight function w*(z) and deduced by taking into account the orthogo- 
nality relation ([8]). 

The secular equation leading to the neutral values of the Rayleigh number can be 
deduced in a similar way by using shifted Legendre polynomials. 
Let 

flo (0, 1) = {/I/, /' G L 2 (0, 1), /(0) = f(l) = 0}, 

be a Hilbert space and denote by the Legendre polynomials defined on (—1,1). 
Then the shifted Legendre polynomials Q (SLP) on (0, 1) are defined by the relation 

Qk( x ) = Lf t (2x — 1) and they are orthogonal on the interval (0, 1), i.e. / QiQjdz = 

Jo 

^ . — - 5ij . Using the identity [9] 

2(2i + l)Qi(z) = Q' i+1 (z) - QUCz), (15) 
we define the set M 2 of orthogonal functions fa, 

fa(z)= [ z Q i (t)dt = Qi+ i7®;r 1 ,i = i,2,... 

Jo 2(2z + 1) 
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that satisfy boundary conditions of the type (f>i(0) = 4>i(l) = at z = and 1 such 
that the set M<i is complete in i?o(0, 1). 

Therefore we can write the unknown functions W, ^, O as series in the form 



•i (z . 



(16) 



i=l 



i=l 



i=l 



The secular equation is obtained following the same steps in the analytical study as 
before, i.e. 



{(D 2 - a 2 )^^) -1 

((D 2 -a 2 )&Ak) -i?a 2 ((l + e/i(z))^,^) 

{(D 2 - a 2 )^,^) {(D 2 -a 2 )^,0 fc ) 



0. 



(17) 



3 Numerical evaluations 

The numerical evaluations were obtained for different significant values of the scale 
parameter e and the wavenumber a. The number of functions in the expansion sets 
was small (n = 4), but the obtained approximative values of the Rayleigh number 
were similar to ones obtained with other methods and we considered that a small 
improvement obtained for n > 4 would not justify more time. In order to compare 
our results and implicitly to test the method, we took into consideration three variable 
gravity fields from [12], i.e. h(z) = —z, h(z) = —z 2 and h(z) = z 2 — 2z. The numerical 
results presented in Tables 1,2,3 show that a decreasing gravity field enlarge the domain 
of stability. For e = 0, we obtained similar evaluations with the classical ones from [2]. 



e 


a 2 


R 2 - SCP 


R 2 - SLP 


0.0 


4.92 


657.512 


675.05 


0.01 


4.92 


660.747 


678.45 


0.03 


4.92 


667.653 


685.33 


0.33 


4.92 


787.363 


808.303 


0.2 


5.00 


730.459 


749.95 


0.2 


9.00 


829.44 


846.70 


0.5 


7.5 


930.982 


952.07 


0.5 


9.00 


994.393 


1015.27 


0.75 


10.0 


1251.178 


1276.05 



Table 1. Numerical values of the Rayleigh number for various values of the parameters for 

h(z) = —z. 
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a 2 


/ i — . ' v / 


R 2 ?rp 

/ i — O ±j± 


n n 


A Q9 




U 1 O.UtJ 


n m 


A Q9 


fi^Q 41 


fi7fi QQ 


u.uo 


A Q9 


UUO.-L / 


fisn on 

UOU.c/U 


u.oo 


A Q9 


795 06 


745 91 


0.2 


5.00 


696.80 


715.87 


0.2 


9.00 


791.24 


808.22 


0.5 


7.5 


813.28 


833.21 


0.5 


9.00 


868.72 


888.54 


0.75 


10.0 


993.51 


1016.20 



Table 2. Numerical values of the Rayleigh number for various values of the parameters for 

h{z) = -z 2 . 



e 


a 2 


K 2 - SCP 


R 2 - SLP 


0.0 


4.92 


657.512 


675.05 


0.01 


4.92 


662.29 


679.91 


0.03 


4.92 


671.95 


689.83 


0.33 


4.92 


861.25 


882.05 


0.2 


5.00 


767.40 


787.44 


0.2 


9.00 


871.37 


889.03 


0.5 


7.5 


1088.2 


1110.46 


0.5 


9.00 


1162.4 


1184.11 


0.75 


10.0 


1687.8 


1713.45 



Table 3. Numerical values of the Rayleigh number for various values of the parameters for 

h(z) = z 2 - 2z. 

4 Conclusions 

In this paper we presented two shifted polynomials-based methods for the study 
of the linear stability of the mechanical equilibrium of a horizontal layer of a viscous 
incompressible fluid heated from below in the case of a variable gravity field. 

We provided numerical results for various gravity fields, decreasing but not all 
linear. These results proved to agree quite well with the results obtained by us with 
the Galerkin method based on trigonometric functions [3] and also with the classical 
one existing in the literature. The numerical evaluations of the Rayleigh number were 
obtained for different values of the physical parameters, allowing a conclusion on the 
effects of these parameters on the stability domain. It was seen that the domain of 
stability decreases as the gravity field is increasing. 

Although both expansion sets of functions led to good numerical evaluations of 
the Rayleigh number, it seems like the method based on the expansion upon shifted 
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Chebyshev polynomials was more effective. However, in most cases, the Legendre 
polynomials are preferred in the Galerkin approach and the Chebyshev polynomials 
are considered suitable for the collocation methods. 
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